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1 Abstract 

Effects of the bond dilution on the critical temperatures, phase diagrams and the magnetization 
behaviors of the isotropic and anisotropic quantum Heisenberg model have been investigated in 
detail. For the isotropic case, bond percolation threshold values have been determined for several 
numbers of two (2D) and three (3D) dimensional lattices. In order to investigate the effect of the 
anisotropy in the exchange interaction on the results obtained for the isotropic model, a detailed 
investigation has been made on a honeycomb lattice. Some interesting results, such as second order 
reentrant phenomena in the phase diagrams have been found. Keywords: Quantum anisotropic 
Heisenberg model; bond dilution; bond percolation threshold 

2 Introduction 

Quenched randomness effects are very important in modeling real materials, since consideration of 
these effects simulates a realistic model of real materials. Real materials have some uncontrollable 
defects and these defects can be modeled by introducing site dilution (randomly distributed non 
magnetic atoms), bond dilution (randomly broken bonds between the magnetic atoms) or both 
of them into the related model. It is a well known fact that, Heisenberg model produces more 
realistic results than the Ising model, in order to explain the magnetic properties of real materials. 
Thus, it is important to work on the Heisenberg model with these quenched randomness effects. 
These quenched randomness effects produce different behaviors in the magnetic properties of the 
model, e.g. different phase transition characteristics from the pure model (i.e. the model without 
any quenched randomness effects). 

Heisenberg model with quenched randomness effects has been studied widely by a variety of 
methods, such as spin-1/2 (S-l/2) isotropic Heisenberg model with bond dilution on 2D lattices 
with Monte Carlo (MC) simulation [1] [2] [3] [4] , anisotropic model (by means of the XXZ model) 
on 2D lattices with real space renormalization group (RSRG) technique [5]. On the other hand, 
random bond distributed systems in which spin glass phases originate have been studied, e.g. 
discrete distribution on S-l/2 Heisenberg model with pair approximation [6] [7], density matrix 
product approximation [5], Gaussian distribution with imaginary time Grassmann field theory [5] 
and S-l Heisenberg model with discretely distributed random bonds with exact diagonalization 
method [TU], Besides, both site and bond diluted systems have been studied with quantum MC 
on 2D lattices [TT] and S-l/2 Heisenberg model with site-bond correlated dilution (which covers 
uncorrelated site dilution as a limit) on 2D and 3D lattices with RSRG [T2] and MC [T3] , and also 
using a variational principle for the free energy [T3]. All of these work are related to the isotropic 
Heisenberg model or XXZ model. 

The aim of this work is determine the effect of the bond dilution on the phase diagrams and the 
thermodynamic properties of the anisotropic quantum Heisenberg model. By anisotropy, we do 
not restrict ourselves with XXZ model. Namely, we want to determine the effect of the anisotropy 
in the exchange interaction on the phase transition characteristics not only by means of XXZ 
model. The method is effective field theory (EFT) with two spin cluster approximation [TB]. EFT 
approximation can provide results that are superior to those obtained within the traditional mean 
field approximation (MFA), due to the consideration of self spin correlations which are omitted in 
the MFA. 

EFT for a typical magnetic system starts by constructing a finite cluster of spins which repre- 
sents the system. Callen-Suzuki spin identities [T7] are the starting point of the EFT for the 
one spin clusters. If one expands these identities with differential operator technique [18], multi 
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spin correlations appear, and in order to avoid from the mathematical difficulties, these multi spin 
correlations are often neglected by using decoupling approximation (DA) |19| . Working with larger 
finite clusters will give more accurate results. Callen-Suzuki identities have been generalized to two 
spin clusters in Ref. [20] (namely EFT-2 formulation). This EFT-2 formulation has been success- 
fully applied to a variety of systems, such as quantum S-l/2 Heisenberg ferromagnet [21] [22] and 
antifcrromagnet [23] systems, classical n- vector model [MIES], and spin-1 Heisenberg ferromagnet 

[23 

The paper is organized as follows: In Sec. [31 we briefly present the model and formulation. 
The results and discussions are presented in Sec. [4] and finally Sec. [5] contains our conclusions. 



3 Model and Formulation 



We consider a lattice which consists of N identical spins (S-l/2) such that each of the spins has z 
nearest neighbors. The Hamiltonian of the system is given by 
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Jy and J\ LJ> stand for the 



where sf, sf and sf denote the Pauli spin operators at a site 
components of the exchange interaction J (in other words anisotropy in the exchange interaction) 
between the nearest neighbor spins i and j . The sum is carried over the nearest neighbors of the 
lattice. The bonds between the spins i and j are randomly distributed via 



P = cS (jM - j) + (1 - c) S (j&A 



(2) 

in the bond dilution problem. The distribution given by Eq. ([2]) distribute bonds randomly between 
lattice sites as c percentage of bonds are closed and remaining 1 — c percentage of bonds are open, 
i.e. c is the concentration of closed bonds in the lattice. Here 6 stands for the delta function and 
c is a real number which is defined within the range of < c < 1. The distribution given by Eq. 
([2"]) reduces to the system with homogenously distributed bonds (i.e. pure system) for c = 1. 

We use the two spin cluster approximation as an EFT formulation, namely EFT-2 formulation [15] 
In this approximation, we choose two spins (namely s\ and S2) and treat the interactions in this 
two spin cluster exactly. In order to avoid some mathematical difficulties we replace the perimeter 
spins of the two spin cluster by Ising spins (axial approximation) [25] . With the procedure defined 
in Ref. [25j . we get an expression for the magnetization per spin as 
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where /? = l/(fcsT), fee is Boltzmann constant and T is the temperature. The inner average 
bracket in Eq. ([3]) (which has no subscript) stands for thermal average and the outer one (which 
has subscript r) is for the configurational averaging which is necessary for including the effect of 
the random bond distribution. The parameters in Eq. ^ are given by 
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where stands for the z component of the nearest neighbor of the spin s\ while sf stands for 
the z component of the nearest neighbor of the spin S2- The sums in the Eq. ^ are over the 
nearest neighbor sites of the sites labeled 1 and 2, respectively. The configurational averages can 
be calculated via integration of the expression by using Eq. ([2]) , over the all bonds of the treated 
cluster. 



dJ (12)p Ul2)\ JJdjWrfjWp (j( lfe >) P (j 

k,l 



h' 

(21) \ n + 



sinh (PX' Q ) 



X o cosh {pX' Q ) + cxp (-2/3 Ji 12) ) cosh {fiY£) 

(7) 



where product is taken over the nearest neighbors of the sites 1 and 2. If we perform the integration 
concerning the bond between the sites 1 and 2 then we get the following expression, 
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Eq. (jSJ) can be written in terms of the differential operators using differential operator technique 
[TBI and it is found as 
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where 
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In Eq. (|10[) . the parameters = d/dx and V B = <9/cfy are the usual differential operators in 
the differential operator technique. Differential operators act on an arbitrary function via 



exp (aV x + bV y )g (x, y) = g (x + a,y + b) 



(13) 



with any constant a and b. 

Now let us assume that each of si and S2 has number of zo distinct nearest neighbors and 
both of them have z± common nearest neighbors. This means that, in Eq. ([B]), the two sums have 
number of Z\ common terms. If we take the integration in Eq. ([TU| with the help of Eqs. © and 
© and using DA [H] we get an expression, 



m = ([A x + mB x ] z ° [A y + mB y } z ° [A xy + mB xy \ Zl ) f (x, y) \ x=0 ^ y=0 
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for the magnetization. The coefficients are defined by 

A x = ccosh(J z V x ) + (1 -c) B x = csmh(J z V x ) 

A y = ccosh(J z V y ) + (1 - c) B y = csmh(J z W y ) (15) 

A xy = ccosh [J z (V x + V y )] + (1 - c) B xy = csinh [J z (V, + V„)]. 

With the help of the Binomial expansion, Eq. ([Til) can be written as 

Zq Zo Z\ 

m = J2J2J2 c ' Pq r^ p+q+r (16) 

p=0<j=0r=0 

where the coefficients are 

C' pqr = ( Z p)( Z q)( "r ) K°- V K~ qA %~ rB l B l B lyf^y) l*=o,v=o (17) 

and these coefficients can be calculated by using the definitions given in Eqs. (JTT|) . ([TB")) and (|T5|) . 
Let us write Eq. (|16p in more familiar form as 

z 

m = Y,C k m k (18) 

k=0 

Zq Zq Z\ 

C k = y^y^y^ ^p+g+r.fcCpqr- (19) 
p=0<j=0r=0 

where Sij is the Kronecker delta. 

For a given set of Hamiltonian parameters (J x , J y , J z ), temperature and bond distribution 
parameter (c), we can determine the coefficients from Eq. (|T71) and we can obtain a non linear 
equation from Eq. (|18[) . By solving this equation, we can get the magnetization (m) for a given 
set of parameters and temperature. Since the magnetization is close to zero in the vicinity of the 
critical point then we can obtain a linear equation by linearizing the equation given in Eq. (|18|) 
which allows us to determine the critical temperature. Since we have not calculated the free energy 
in this approximation, we can locate only second order transitions from the condition given as 

d = 1, C 3 < 0. (20) 

The tricritical point at which second and first order transition lines meet can be determined from 
the condition 

d = l, C 3 =0. (21) 



4 Results and Discussion 

Let us scale the exchange interaction components with the unit of energy J as 

r n = y (22) 

where n = x,y, z. Let us choose r z = 1, then r x , r y can be used as the measure of the anisotropy 
in the exchange interaction. It can be seen from the formulation of the problem (e.g. see Eqs. 
(fTTj) - (fTS"|) ) that the transformation J x —> J y ,J y — > J x does not change the expression of the 
order parameter, regardless of the bond distribution. Hence, we can say that aforementioned 
transformation does not affect the formulation. Therefore, let us fix r x and concentrate only on 
the varying r y values. Our investigation will be focused on selected several 2D and 3D lattices 
in isotropic case and on a honeycomb lattice in anisotropic case with the ferromagnetic exchange 
interaction, i.e. r x , r y , r z > 0. It can be seen from the Hamiltonian of the system which is given by 
Eq. UJ) that for r z > 0, r x = r y = the system reduces to a model which has Ising type symmetry 
and r x = 0, r y , r z > or r x , r z > 0, r y = reduces the system to the system that have XY type 
symmetry. 
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4.1 Isotropic Model 



In this case, all components of the exchange interaction are equal to each other, i.e. r n = 1 in Eq. 
(|22l) . In order to represent the behavior of the critical temperature (/cbT c /J) as a function of the 
bond concentration (c) for the isotropic quantum Heisenberg model, we depict the phase diagrams 
for three different lattices in (fc^Tc/J, c) plane in Fig. (H])(a). As we can see from Fig. (HJ(a) that 
the critical temperature of the system gradually decreases when the bond concentration decreases 
and takes the value of zero at a specific value of bond concentration, namely at the bond percolation 
threshold (c*) which depends on the lattice geometry. At this point, we have to mention that bond 
percolation threshold value may not correspond to the value that makes the critical temperature 
zero. If the phase diagram of the system exhibits a reentrant behavior in the (/cbT c /J, c) plane 
then the bond percolation threshold value corresponds to the c value that terminates the phase 
diagram, i.e. for the value that provide c < c* system can not reach the ordered phase any more. In 
the isotropic case of the quantum Heisenberg model, the whole phase diagrams in the (fcsT c / J, c) 
plane do not exhibit reentrant behavior. 

The numerical values of bond percolation threshold values obtained within the present work 
are c* = 0.5903 for a honeycomb lattice, c* = 0.4427 for a square lattice, and c* = 0.2974 for a 
simple cubic lattice. We can also see from Fig. (H])(a) that at a certain c, simple cubic lattice has 
highest critical temperature value which is due to the large number of the nearest neighbor spins 
in the lattice. Besides, in Fig. ©(b) we can see the variation of the ground state magnetization 
(mo) with the bond concentration for three lattice geometries. These values of magnetizations 
calculated at the temperature fc^T/J = 0.001 can be treated as ground states since the energy 
comes from the thermal fluctuations are negligible in comparison with the energy comes from the 
spin-spin interaction. As we can see from Fig. (JTJ) (b) that, ground state magnetization starts from 
the value of too — 1.0 and, stays in this value for a while when the bond concentration decreases. 
After then, it starts to fall and reach the zero at a bond percolation threshold value. 




Figure 1: (a) The variation of the critical temperature (ksTc/J) of the system with the bond 
concentration (c), (b) the variation of the ground state magnetization of the system with the bond 
concentration (c), for selected lattices with the isotropic quantum Heisenberg model. The fixed 
parameter values are r x = r y = r z = 1.0. In (b), the temperature is selected as ksTc/J = 0.001. 

It is a known fact that for a certain set of Hamiltonian parameters the critical temperature 
of the pure system with Ising type symmetry is greater than the system with XY type symmetry 
and this value is greater than that of the system with Heisenberg type symmetry [15]. From this 
fact we expect that, the relation between the bond percolation threshold values of three types 
of symmetries to be c*(H) > c*{XY) > c*(I). Here, the letters in the parenthesis depict the 
symmetry type of the Hamiltonian. Namely, H stands for Heisenberg type symmetry (with the 
anisotropy values r x — r y = r z = 1.0), XY stands for XY type symmetry (with the anisotropy 
values r x = 0.0, r y = r z = 1.0), and / stands for Ising type symmetry (with the anisotropy values 
r x = r y = 0.0, r z = 1.0). But the formulation gives this relation as c*(XY) > c*(H) > c*(I), as 
seen in Table (1). We will touch upon this point in the next subsection. The bond percolation 
threshold values of the different lattices for different Hamiltonian symmetries can be seen in Table 
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(1). 



Tabic 1: The bond percolation threshold values for several 
lattices with the Hamiltonians that have Ising type symmetry 
(I), XY type symmetry (XY) and Hcisenberg type symmetry 
(H). 
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The effect of the bond concentration on the behavior of the magnetization with temperature 
can be seen in Figs, ©(a) for a honeycomb lattice and ©(b) for a simple cubic lattice. As seen 
in Fig. ([2]) that, decreasing c causes a decline in the critical temperature value, also after some 
value of it, it causes decline in the ground state magnetization. Ground state magnetization starts 
to decrease at a lower value of c for the simple cubic lattice (in comparison with the honeycomb 
lattice), as seen in Fig. ([2]) in compatible with Fig. (H])(b) (e.g. all curves of c value that provide 
0.6 < c < 1.0 for the simple cubic lattice have ground state magnetization ?n = 1.0 while for 
the honeycomb lattice the same ground state magnetization value has been observed only for 
c = 1.0, 0.9 values). Again, this is due to the excess number of nearest neighbor spins of the simple 
cubic lattice in comparison with those of the honeycomb lattice. 
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Figure 2: Variation of magnetization with temperature for (a) honeycomb and (b) simple cubic 
lattices for selected bond concentration values within the quantum isotropic Heisenberg model. 
The fixed parameter values are r x = r y = r z = 1.0. 



4.2 Anisotropic Model 

Now, let us see the effect of the anisotropy in the exchange interaction on the results obtained 
for the isotropic case in the preceding subsection. To do this, we start with the phase diagrams 
corresponding Fig. ([T]) (a) for the isotropic case. The phase diagrams for a honeycomb lattice in 
(fcsT c /J, c) plane can be seen in Fig. §5$ for different anisotropy values. In Fig. ©(a), we can 
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Figure 3: Variation of the critical temperature (ksTc/ J) of the anisotropic quantum Heisenberg 
model as a function of bond concentration (c) for a honeycomb lattice, with some selected values 
of anisotropy in the exchange interaction (r x ,r y ). The fixed parameter value is r z = 1.0 



see the phase diagrams of the system with Ising type symmetry (the curve labeled by A) and XY 
type symmetry (curves labeled by B,C,D). As we can see from the curves labeled by B,C,D in Fig. 
([3]) (a), as the amount of anisotropy increases then the ferromagnetic region gets narrower in the 
(kBT c /J, c) plane for the XY type model and after a certain value of c, phase diagrams disappears. 
This value depends on the lattice geometry. Bond percolation threshold values of the XY type 
system gradually decrease when the anisotropy increases. There is not any reentrant behavior 
observed in the phase diagrams of the XY type system. In Figs. (J3|) (b) and (c), transition from 
the curves labeled by A to the other curves corresponds to the transition from the system with 
XY type symmetry to the system with H type symmetry. When we take a look at Fig. ([3]) (b), 
we see that increasing anisotropy effects the phase diagrams in the (fcsT c /J, c) plane in a more 
complicated way in comparison with Fig. ([3]) (a). Firstly, as the anisotropy increases then the 
bond percolation thresholds decrease for a while then tend to increase. We can clearly see here, 
although the critical temperature of the pure system for the XY type symmetry (the curve labeled 
A in Fig. ©(b)) is greater than that of the H type symmetry (curve labeled C in Fig. (0(b)), 
this docs not imply that the bond percolation threshold value of the XY type symmetric model to 
be lower than the H type symmetric model. This behavior of the curves labeled A and C seen in 
Fig. (|3])(b) verifies the results in Table (1). Secondly, increasing anisotropy can induce a reentrant 
behavior in the phase diagrams (e.g. see the curve labeled by B in Fig. © (b)). Besides, we can 
see from Fig. ([3]) (c) that, increasing anisotropy in the exchange interaction can produce phase 
diagrams that behave qualitatively different from the others. For instance, for the curve labeled 
by E in Fig. ([3]) (c), decreasing anisotropy first causes a sharp decrease in the critical temperature 
then the critical temperature stays almost constant before it depresses to zero. In this figure, the 
curves labeled by B and D also exhibit reentrant behavior. One other interesting behavior can be 
seen in Fig. ^ (d). For r x = 3.0, with low anisotropics (i.e. small values of r y ), there arc no 
phase diagrams. The system can not exhibit an ordered phase. When the value of r v increases 
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then the phase diagrams appear and then disappear again. In conclusion, we can say that while the 
transition from the I type symmetry to the XY type symmetry (Fig. (a)) rising anisotropy in 
the exchange interaction can not induce reentrant behavior in the phase diagrams which are in the 
(ksTc/ J, c) plane. But, if increasing anisotropy causes a transition from the XY type symmetry 
to the H type symmetry then it can induce a reentrant behavior in the phase diagrams in the 
(k B T c /J,c) plane, (Fig. © (b) and (c)). 

In order to concentrate on the evolution of the bond percolation threshold value with anisotropy 
in the exchange interaction, let us investigate the variation of c* value with r y for some selected 
values of r x . These curves can be seen in Fig. (jl])(a). The curve labeled by A corresponds to 
the XY type system. As mentioned in the comments of Fig. Q, for the XY type system, with 
increasing anisotropy, the bond percolation threshold value smoothly increases to the value 1.0 
while for the system with H symmetry, the bond percolation threshold value first decreases with 
increasing anisotropy then increases and reaches to the value of 1.0 (see the curves labeled by 
B,C,D in Fig. 0(a)). When r x increases, the c* value corresponding to r y = increases and 
after a specific value of r x = 2.7105 the c* value for r y = becomes 1.0. This value is nothing 
but the value of the curve labeled by A intersects the c* = 1.0, due to the symmetry of the 
formulation as mentioned at the beginning of this section. This is just the value that makes the 
critical temperature of the pure system (i.e. c = 1.0) equal to zero with XY type symmetry. We 
can numerically determine the similar borders of the model with H type symmetry. 

As seen in Fig. Q (a), curves intersect the c* = 1.0 line at a specific value of r y which depends 
on the value of r x . The curves labeled A,B and C intersect the c* = 1.0 line only by one point while, 
the curve labeled D by two point. Let denote these two intersection points as r* and r* which 
provide r* > r* . In other words r * denotes the value of r y which is related r x curve intersects to 
the c* = 1 line on the right side and r* denotes the value of r y which related r x curve intersects to 
the c* = 1 line on the left side, if present. These two special values in the (r x , r y ) plane can be seen 
in Fig. (j4]) (b). For any point which is outside of the region confined by this curve and r x , r y axes, 
we can say that the system can not exhibit an ordered phase for a honeycomb lattice. This border 
consists of three segments: r y = r x + 2.7105 line, r y = r x — 2.7105 line and the remaining curve. 
The region that the system can exhibit an ordered phase becomes expanded when r x increases 
from 0.0 to 1.0, after then that region starts to shrink. We can say from the symmetry properties 
of the formulation that, this border curve is symmetric about the r y = r x line. We can obtain 
similar borders for the system which has bond concentration values different from c = 1.0. These 
curves can be seen in Fig. (j4|) (c). Again we can say with Fig. ((4]) (c) that, for any point which 
is outside of the region confined by the related c curve and r x ,r y axes, the system with (1 — c) 
percentage open bonds, can not exhibit an ordered phase for a honeycomb lattice. 

Finally, we want to investigate the effect of the bond dilution on the ground state magnetizations 
of the system. In the isotropic case, this effect was as decreasing c gradually decreases the ground 
state magnetizations (Fig. (p} (b)). But in anisotropic case, this situation is a little complicated 
as seen in Fig. ([5]), which is the variation of the ground state magnetization (mo) with anisotropy 
(r y ) and bond concentration (c) for some selected values of r x . Fig. (|S|) (a) is just concerning the 
model with XY type symmetry. There is not any complicated behavior as in the isotropic model. 
For a certain c, the ground state magnetization regularly decreases with increasing r y . 

But the situation is more complicated in Figs. ([5]) (b) and (c). These are the surfaces regarding 
the model with H type symmetry and all the reentrant behaviors in the phase diagrams of the 
system (e.g. curve labeled by B in Fig. © (b) and (c)) shows itself in the behavior of the ground 
state magnetizations. Beside these, r x = 2.0 surface is not smooth while r x = 1.0 surface is 
smooth i.e., for r x = 1.0 with any value of c, increasing r y changes the ground state magnetization 
continuously but this is not the case for r x = 2.0. Increasing r x value also shrinks this surface as 
expected. Fig. ©(d) has different characteristics than the others; the ground state magnetization 
surface can not reach the small valued r y regions, as well as the higher valued r y regions. This fact 
is compatible with the curve labeled by D in Fig. (|4|) (a). 



X 




Figure 4: (a) Variation of the bond percolation threshold value (c*) with r y for some selected 
values of r x , (b) the border of the ferromagnetic region of the pure (c = 1.0) model in the (r x , r y ) 
plane, (c) equally valued bond percolation threshold curves in the (r x ,r y ) plane for a honeycomb 
lattice within anisotropic quantum Hciscnbcrg model. The fixed parameter value is r z = 1.0. 

5 Conclusion 

The effect of the bond dilution on the phase diagrams of the isotropic and anisotropic quantum 
Heisenberg model has been investigated in detail. First in the isotropic model, effect of the bond 
dilution process on the critical temperatures have been given for honeycomb, square and simple 
cubic lattices. There is not any reentrant behavior observed in that phase diagrams plotted in 
(fesTc/J, c) plane for the isotropic Heisenberg model. Also, the variation of the ground state 
magnetization with bond concentration (c) and the evolution of the variation of the magnetization 
with temperature curves with c have been presented. The bond percolation threshold values which 
is the value of c that terminates the phase diagrams in the (kBT c / J,c) plane have been given for 
several 2D and 3D lattices with the Hamiltonians that have different symmetries. 

In addition, the same investigation has been performed for the anisotropic model. Anisotropic 
model has been treated within three different regions. Namely, the Hamiltonian that have Ising type 
symmetry r x = r y = 0.0, r z = 1.0, XY type symmetry r x = 0.0, r y > 0.0, r z = 1.0 and Heisenberg 
type symmetry r x ,r y > 0.0, r z = 1.0. All these models that have different symmetry properties 
have been investigated in detail: The effect of the bond dilution process on both the ground state 
magnetizations and critical temperatures have been determined via the phase diagrams in the 
{UbTc/ J,c) plane. Beside these, the bond percolation threshold values (c*) have been determined 
for the anisotropic model and depicted in the (c*,r y ) plane for selected values of r x . For the XY 
type symmetric model, decreasing bond concentration results in gradually increase of the bond 
percolation threshold value. There is not any reentrant behavior observed within the model that 
have XY type symmetry on the phase diagrams in the (/csT c /J, c) plane. 

On the other hand, in the model that have Hciscnbcrg type symmetry, this situation is more 
complicated. Increasing anisotropy does not results in a gradually increase of the c* value as in 
the model that have XY type symmetry. However, increasing anisotropy first causes a decrease in 
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Figure 5: Variation of the ground state magnetization with c and r v for a honeycomb lattice 
within the quantum Hciscnbcrg model for some selected values of r x . The fixed parameter value 
is r x = 1.0. 

the c* value then an increase in it before it reaches the value of 1.0. Beside this fact, increasing 
anisotropy also give rise to the reentrant behavior in the phase diagrams in the (fcsT c /J, c) plane. 
All of these reentrant behaviors are found to be of the second order. 

We hope that the results obtained in this work may be beneficial form both theoretical and 
experimental point of view. 
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